struct Vector {
	double x, y, z;

	Vector (): x(0), y(0), z(0) {}

	Vector (double _x, double _y, double _z): x(_x), y(_y), z(_z) {}

	Vector (float* feedbackBuffer): 
		x	( feedbackBuffer[1] ), 
		y	( feedbackBuffer[2] ), 
		z	( feedbackBuffer[3] ) 
	{
//		cout << toString().c_str() << endl;	
	}

	~Vector () {}

	bool	operator == (const Vector& v) const { return x == v.x && y == v.y && z == v.z;}

	Vector	operator -	()				  const	{ return Vector (-x, -y, -z); } 
	Vector	operator -	(const Vector& v) const	{ return Vector (x - v.x, y - v.y, z - v.z); }
    void	operator -= (const Vector& v)		{ x -= v.x; y -= v.y; z -= v.z; }
	void	operator -= (double d)				{ x -= d; y -= d; z -= d; }

	Vector	operator +	(const Vector& v) const { return Vector(x + v.x, y + v.y, z + v.z); }
	void	operator += (const Vector& v)		{ x += v.x; y += v.y; z += v.z; }
	void	operator += (double d)				{ x += d; y += d; z += d; }

	Vector	operator *	(const double& d) const { return Vector(x*d, y*d, z*d); }
	void	operator *= (const double& d)		{ x *= d; y *= d; z *= d; }

	Vector	operator /	(double d)		  const { d = 1/d; return Vector(x*d, y*d, z*d); }
	void	operator /= (double d)				{ d = 1/d; x *= d; y *= d; z *= d; }

	double	length			() const { return sqrt( x*x + y*y + z*z ); }
	double	lengthSquared	() const { return x*x + y*y + z*z; }
	void	normalize		()		 { if (x != 0 || y != 0 || z != 0) *this /= length(); }

	double	dot		(const Vector& v) const { return x*v.x + y*v.y + z*v.z; }
	double	angle	(const Vector& v) const { return acos( dot(v) / (length() * v.length()) ); }

	Vector	cross	(const Vector& v) const { 
		return Vector( (y * v.z) - (z * v.y), 
					   (z * v.x) - (x * v.z), 
					   (x * v.y) - (y * v.x) ); 
	}

	string toString() const {
		char buffer[255];
		sprintf(buffer, "%+f %+f %+f", x, y, z);
		return string(buffer);
	}
	//

};
//

typedef Vector Point;
typedef	Vector Vertex;
//


struct Plane {
	Vector normal;
	double k;

	bool vertOutside (Vector vert) {
		if (vert.dot(normal) < k)
			return false;

		return true;
	}
	//

};



class Transformation {
public:
	// these are stored in column-major,  for OpenGL
	// compatibility, but operations are expressed in 
	// row-major, 'cause that's how I think. 
	double cell[16];

	Transformation () {
		setToIdentity();
	}
	//

	Transformation(
		double m11, double m12, double m13, double m14,
		double m21, double m22, double m23, double m24,
		double m31, double m32, double m33, double m34,
		double m41, double m42, double m43, double m44) {
			
		set(m11, m12, m13, m14, 
		    m21, m22, m23, m24, 
		    m31, m32, m33, m34, 
		    m41, m42, m43, m44);
	}
	//

	~Transformation () {}

	void set (
		double m11, double m12, double m13, double m14, 
		double m21, double m22, double m23, double m24, 
		double m31, double m32, double m33, double m34,
		double m41, double m42, double m43, double m44) {
		cell[0] = m11; cell[4] = m12; cell[8]  = m13; cell[12] = m14;
		cell[1] = m21; cell[5] = m22; cell[9]  = m23; cell[13] = m24;
		cell[2] = m31; cell[6] = m32; cell[10] = m33; cell[14] = m34;
		cell[3] = m41; cell[7] = m42; cell[11] = m43; cell[15] = m44;
	}
	//

	void setToIdentity () {
		set( 1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 );
	}

	void loadFromStack (int mode = currentMatrixStack()) {
		glGetDoublev (mode, cell);	
	}

	void loadToStack (GLenum mode) {
		int savedMode; glGetIntegerv (GL_MATRIX_MODE, &savedMode);
		glMatrixMode  (mode);
		glLoadMatrixd (cell);
		glMatrixMode  (savedMode);
	}
	//

	void transpose () {
		set( cell[0],  cell[1],  cell[2],  cell[3], 
			 cell[4],  cell[5],  cell[6],  cell[7],
			 cell[8],  cell[9],  cell[10], cell[11],
			 cell[12], cell[13], cell[14], cell[15] );
	}
	//

	Vertex operator * (const Vertex& v) const {
		double x = v.x * cell[0]  +  v.y * cell[4]  +  v.z * cell[8]   +  cell[12];
		double y = v.x * cell[1]  +  v.y * cell[5]  +  v.z * cell[9]   +  cell[13];
		double z = v.x * cell[2]  +  v.y * cell[6]  +  v.z * cell[10]  +  cell[14];
		double w = v.x * cell[3]  +  v.y * cell[7]  +  v.z * cell[11]  +  cell[15];

		return Vertex (x, y, z)	/ w  ;
	}
	//

	Transformation operator * (const Transformation& a) const {
		return Transformation (
			cell[0] * a.cell[0]   +  cell[4] * a.cell[1]   +  cell[8] * a.cell[2]   +  cell[12] * a.cell[3],
			cell[0] * a.cell[4]   +  cell[4] * a.cell[5]   +  cell[8] * a.cell[6]   +  cell[12] * a.cell[7],
			cell[0] * a.cell[8]   +  cell[4] * a.cell[9]   +  cell[8] * a.cell[10]  +  cell[12] * a.cell[11],
			cell[0] * a.cell[12]  +  cell[4] * a.cell[13]  +  cell[8] * a.cell[14]  +  cell[12] * a.cell[15],

			cell[1] * a.cell[0]   +  cell[5] * a.cell[1]   +  cell[9] * a.cell[2]   +  cell[13] * a.cell[3],
			cell[1] * a.cell[4]   +  cell[5] * a.cell[5]   +  cell[9] * a.cell[6]   +  cell[13] * a.cell[7],
			cell[1] * a.cell[8]   +  cell[5] * a.cell[9]   +  cell[9] * a.cell[10]  +  cell[13] * a.cell[11],
			cell[1] * a.cell[12]  +  cell[5] * a.cell[13]  +  cell[9] * a.cell[14]  +  cell[13] * a.cell[15],

			cell[2] * a.cell[0]   +  cell[6] * a.cell[1]   +  cell[10] * a.cell[2]   +  cell[14] * a.cell[3],
			cell[2] * a.cell[4]   +  cell[6] * a.cell[5]   +  cell[10] * a.cell[6]   +  cell[14] * a.cell[7],
			cell[2] * a.cell[8]   +  cell[6] * a.cell[9]   +  cell[10] * a.cell[10]  +  cell[14] * a.cell[11],
			cell[2] * a.cell[12]  +  cell[6] * a.cell[13]  +  cell[10] * a.cell[14]  +  cell[14] * a.cell[15],

			cell[3] * a.cell[0]   +  cell[7] * a.cell[1]   +  cell[11] * a.cell[2]   +  cell[15] * a.cell[3],
			cell[3] * a.cell[4]   +  cell[7] * a.cell[5]   +  cell[11] * a.cell[6]   +  cell[15] * a.cell[7],
			cell[3] * a.cell[8]   +  cell[7] * a.cell[9]   +  cell[11] * a.cell[10]  +  cell[15] * a.cell[11],
			cell[3] * a.cell[12]  +  cell[7] * a.cell[13]  +  cell[11] * a.cell[14]  +  cell[15] * a.cell[15]);
	}
	//

	void operator *= (const Transformation& a) {
		//this = this * a
		double r11 = cell[0] * a.cell[0]   +  cell[4] * a.cell[1]   +  cell[8] * a.cell[2]   +  cell[12] * a.cell[3];
		double r12 = cell[0] * a.cell[4]   +  cell[4] * a.cell[5]   +  cell[8] * a.cell[6]   +  cell[12] * a.cell[7];
		double r13 = cell[0] * a.cell[8]   +  cell[4] * a.cell[9]   +  cell[8] * a.cell[10]  +  cell[12] * a.cell[11];
		double r14 = cell[0] * a.cell[12]  +  cell[4] * a.cell[13]  +  cell[8] * a.cell[14]  +  cell[12] * a.cell[15];

		double r21 = cell[1] * a.cell[0]   +  cell[5] * a.cell[1]   +  cell[9] * a.cell[2]   +  cell[13] * a.cell[3];
		double r22 = cell[1] * a.cell[4]   +  cell[5] * a.cell[5]   +  cell[9] * a.cell[6]   +  cell[13] * a.cell[7];
		double r23 = cell[1] * a.cell[8]   +  cell[5] * a.cell[9]   +  cell[9] * a.cell[10]  +  cell[13] * a.cell[11];
		double r24 = cell[1] * a.cell[12]  +  cell[5] * a.cell[13]  +  cell[9] * a.cell[14]  +  cell[13] * a.cell[15];

		double r31 = cell[2] * a.cell[0]   +  cell[6] * a.cell[1]   +  cell[10] * a.cell[2]   +  cell[14] * a.cell[3];
		double r32 = cell[2] * a.cell[4]   +  cell[6] * a.cell[5]   +  cell[10] * a.cell[6]   +  cell[14] * a.cell[7];
		double r33 = cell[2] * a.cell[8]   +  cell[6] * a.cell[9]   +  cell[10] * a.cell[10]  +  cell[14] * a.cell[11];
		double r34 = cell[2] * a.cell[12]  +  cell[6] * a.cell[13]  +  cell[10] * a.cell[14]  +  cell[14] * a.cell[15];

		double r41 = cell[3] * a.cell[0]   +  cell[7] * a.cell[1]   +  cell[11] * a.cell[2]   +  cell[15] * a.cell[3];
		double r42 = cell[3] * a.cell[4]   +  cell[7] * a.cell[5]   +  cell[11] * a.cell[6]   +  cell[15] * a.cell[7];
		double r43 = cell[3] * a.cell[8]   +  cell[7] * a.cell[9]   +  cell[11] * a.cell[10]  +  cell[15] * a.cell[11];
		double r44 = cell[3] * a.cell[12]  +  cell[7] * a.cell[13]  +  cell[11] * a.cell[14]  +  cell[15] * a.cell[15];
		set (r11,r12,r13,r14,r21,r22,r23,r24,r31,r32,r33,r34,r41,r42,r43,r44);
	}
	//

	void setToTranslation(double x, double y, double z) {
		set ( 1,0,0,x, 0,1,0,y, 0,0,1,y, 0,0,0,1 );
	}

	void setToTranslation(const Vertex& v) {
		setToTranslation(v.x, v.y, v.z);
	}

	void setToRotation(double angle, Vertex axis){
		axis.normalize();
		angle *= 3.14159265 / 180;
		double cosine = cos(angle);
		double sine   = sin(angle);

		cell[0] = axis.x * axis.x * (1 - cosine) + cosine;
		cell[4] = axis.x * axis.y * (1 - cosine) - (axis.z * sine); 
		cell[8] = axis.x * axis.z * (1 - cosine) + (axis.y * sine); 
		cell[12] = 0;

		cell[1] = axis.y * axis.x * (1 - cosine) + (axis.z * sine); 
		cell[5] = axis.y * axis.y * (1 - cosine) + cosine;  
		cell[9] = axis.y * axis.z * (1 - cosine) - (axis.x * sine);  
		cell[13] = 0;

		cell[2] = axis.z * axis.x * (1 - cosine) - (axis.y * sine); 
		cell[6] = axis.z * axis.y * (1 - cosine) + (axis.x * sine); 
		cell[10] = axis.z * axis.z * (1 - cosine) + cosine; 
		cell[14] = 0;

		cell[3] = 0; cell[7] = 0; cell[11] = 0; cell[15] = 1;
	}
	//

	void setToScale(double x, double y, double z) {	
		set ( x,0,0,0, 0,y,0,0, 0,0,z,0, 0,0,0,1 );	
	}

	void setToScale(const Vertex& v) {	
		setToScale(v.x, v.y, v.z);
	}

	void scale(double x, double y, double z){
		Transformation T;
		T.setToScale(x, y, z);
		*this *= T;	
	}

	void scale(double factor){
		scale(factor, factor, factor);
	}

	void scale (Vector factor) {
		scale(factor.x, factor.y, factor.z);
	}

	void rotate	(Vector rotation) {
		Transformation T;

		T.setToRotation(rotation.x, Vector(1, 0, 0));
		*this *= T;

		T.setToRotation(rotation.y, Vector(0, 1, 0));
		*this *= T;

		T.setToRotation(rotation.z, Vector(0, 0, 1));
		*this *= T;
	}
	//

	void rotate(double x, double y, double z) {
		rotate( Vector(x, y, z) );
	}

	void rotate (double angle, Vector axis) {
		Transformation T;
		T.setToRotation(angle, axis);
		*this *= T;
	}

	void translate (double x, double y, double z) {
		Transformation T;
		T.setToTranslation(x, y, z);
		*this *= T;			
	}

	void translate (Vector translation) {
		translate(translation.x, translation.y, translation.z);
	}

	void print() const {
		printf("\t %+f %+f %+f %+f\n", cell[0], cell[4], cell[8],  cell[12]);
		printf("\t %+f %+f %+f %+f\n", cell[1], cell[5], cell[9],  cell[13]);
		printf("\t %+f %+f %+f %+f\n", cell[2], cell[6], cell[10], cell[14]);
		printf("\t %+f %+f %+f %+f\n", cell[3], cell[7], cell[11], cell[15]);
		cout << endl;
	}
	//

};